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Abstract 

Sector decomposition is a constructive method to isolate divergences from pa- 
rameter integrals occurring in perturbative quantum field theory. We explain the 
general algorithm in detail and review its application to multi-loop Feynman param- 
eter integrals as well as infrared divergent phase-space integrals over real radiation 
matrix elements. 



1 Introduction 



Modern particle physics has reached a level of experimental accuracy in the percent range, 
and some present and future precision experiments require theoretical uncertainties to 
be at the permille level. This need for precise theory predictions has pushed forward the 
frontier for calculations of higher orders in perturbation theory considerably in recent 
years. The calculation of higher order corrections relies to a large extent on tree- or loop 
level Feynman diagrams, where the unobserved degrees of freedom, respectively the loop 
momenta, have to be integrated out. It is well known that these integrations become 
increasingly difficult at higher orders, as the dimensionality of the integration parameter 
space and/or the number of scales involved is growing. The intricacy is closely related 
to the fact that these integrals in general contain ultraviolet (UV) or infrarecG (IR) 
divergences which need to be renormalised respectively factorised, because they hinder 
an immediate numerical evaluation of complicated expressions. The subtraction of these 
singularities becomes more and more cumbersome at higher orders, due to the fact that 
the divergences will be entangled in an increasingly complicated way. 

For ultraviolet divergences, Bogoliubov, Parasiuk, Hepp and Zimmermann [1-3], es- 
tablished a subtraction scheme valid to all orders in perturbation theory. In fact, the 
original idea of sector decomposition goes back to the proof of the BPHZ theorem by 
Hepp [2] , who used a decomposition of integration parameter space into certain sectors 
in order to disentangle overlapping ultraviolet singularities. 

Concerning infrared divergences, the finiteness of sufficiently inclusive observables is 
guaranteed by the Kinoshita-Lee-Nauenberg theorem [4,5], but in most practical applica- 
tions, more exclusive information about the final state is required, such that a subtraction 
scheme for IR divergences has to be established. A constructive scheme to do so to all 
orders in perturbation theory has not been formulated in full generality yet, at least in 
what concerns soft and collinear divergences. However, the sector decomposition algo- 
rithm as outlined below offers a constructive method to isolate these divergences from 
Feynman parameter integrals, which in principle is valid to all orders in perturbation 
theory. 

For the subtraction of (soft) IR divergences, several approaches have been suggested. 
Early ones, which are based already on the subdivision of the integration parameter 
space into different sectors where the parameters go to zero in an ordered way, can be 
found e.g. in [6-10]. These sectors, associated with certain sets of one-particle irreducible 
subgraphs, are more advanced than the sectors needed in the UV case, and provide a 
resolution of the singularities without the need for an iterative procedure for diagrams 
with off-shell external momenta. Further, the so-called -Reoperation [11, 12] has been 
developed, which removes not only UV-divergences but also soft IR divergences by a 
procedure similar to the i?-operation [1] in the UV case. For a review of these methods, 
see e.g. [13, 14]. 

Later, the decomposition into sectors has been employed to extract logarithmic mass 
singularities from massive multi-scale integrals in the high energy limit at two loops [15]. 
In [16], the concept of sector decomposition has been elaborated to a general algorithm 
in the context of dimensional regularisation, allowing to isolate ultraviolet as well as 
infrared singularities from parameter integrals in an automated way [17]. The algorithm 
also has been implemented in a public code available from [18]. 

1 We will use "infrared" to denote both soft and collinear divergences. 
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At next-to-leading order (NLO) in perturbation theory, the use of dimensional regu- 
larisation [19,20] and the universal infrared structure of gauge theories involving massless 
particles, like QED and QCD, allowed to establish a framework to isolate IR singularities 
analytically, leading to a multitude of phenomenological predictions. More recently, an 
increasing number of results at next-to-next-to- leading order (NNLO) has also become 
available, and quite a number of them made use of the sector decomposition technique 
to isolate infrared singularities or to check analytical results [21-42]. 

As most problems encountered in the calculation of perturbative higher order cor- 
rections can be reduced to the evaluation of multi-dimensional parameter integrals, the 
applicability of sector decomposition is quite universal. 

1.1 Multi-loop diagrams 

The first application of the algorithm presented in [16] was the numerical evaluation of 
massless two-loop box diagrams at certain Euclidean pointfl Pioneering the analytic 
calculation of two-loop box diagrams, the planar topology with light-like external legs 
has been calculated by V.A. Smirnov [43], the non-planar one by J.B.Tausk [44], where 
the numerical results obtained by sector decomposition served as an important check of 
the analytical calculation. For massless two-loop box diagrams with one off-shell leg, 
the numerical results of Ref. [16] were predictions which played a major role to validate 
the subsequent analytical calculations [45,46], finally allowing the calculation of the full 
two-loop QCD matrix element for e + e~ — > 3 jets [47]. 

Subsequently, sector decomposition was used to check a considerable number of an- 
alytical results for two-loop [27,34,35,48-52], three-loop [53-55] and four-loop [36,51] 
diagramajl 

Sector decomposition also has been combined with other methods for an efficient 
numerical calculation of one-loop multi-leg amplitudes, first on a diagrammatic level in 
Refs. [56,57], later for whole amplitudes in Refs. [58,59]. The latter approach contains 
a combination of sector decomposition and contour deformation [60-63], which allows 
to integrate the Feynman parameter representation of an amplitude numerically in the 
physical region. It also contains the application of sector decomposition to phase space 
integrals, which will be discussed below. Similar ideas, i.e. the combination of sector 
decomposition and contour deformation, are employed in Refs. [42,64] for the numerical 
evaluation of multi-loop Feynman diagrams with infrared and threshold singularities. 

Ref. [25] describes the implementation of an algorithm based on sector decomposition 
to extract the 1/e poles as well as the large logarithms of type ln(s/M 2 ) in the high- 
energy limit, allowing to obtain the next-to- leading logarithmic electroweak corrections 
of multi-loop diagrams. 

Despite its success in practical applications, until very recently there was no formal 
proof that a strategy for the iterated sector decomposition can always be found such that 
the iteration is guaranteed to stop. This gap has been filled in Ref. [18], by mapping 
the problem to Hironaka's Polyhedra game [65]. The findings of Ref. [18] subsequently 
have been used to prove that the coefficients of the Laurent scries representing Feynman 
integrals in the Euclidean region with rational values for all invariants are a special class 
of numbers known by mathematicians as periods [66] . 

2 By "Euclidean" we mean that all kinematic invariants formed from external momenta are negative. 
s In Refs. [27,34—36] the checks have been performed with an implementation of sector decomposition 
by M. Czakon, independent from the one in [16]. 
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Also on a more formal level, it is outlined in [67] that sector decomposition can be 
used to automate the renormalisation of quantum field theories, by mapping to the so- 
called Henge decomposition [68], which served to provide a simple proof of the BPHZ 
theorem [69]. This mapping also allowed to establish a direct correspondence between 
overlapping divergences in Feynman parameter space and in momentum space [67]. In 
this context, one should also mention the formulation of renormalisation using Hopf 
algebras [70], which provide a framework to describe the disentanglement of divergent 
subgraphs of Feynman diagrams. 



1.2 Phase space integrals 

After the results for various two-loop box diagrams had become available, the bottleneck 
to make progress in the calculation of differential NNLO cross sections for 1 — > 3 or 2 — > 2 
processes was the complicated infrared singularity structure of phase space integrals over 
matrix elements for the real radiation of (doubly) unresolved massless particles. As these 
phase space integrals can be written as dimensionally regularised parameter integrals, 
sector decomposition can serve to factorise entangled singularity structures in the case 
of real radiation as well. This idea has first been presented in [71] and subsequently 
has been applied to calculate all master four-particle phase space integrals where up 
to two particles in the final state can become soft and/or collinear [21]. Shortly after, 
this approach has been extended to be applicable to exclusive final states as well by 
expressing the functions produced by sector decomposition in terms of distributions [22]. 
A rapid development [23,24] lead to the calculation of e + e~ — > 2 jets at 0{a 2 s ) [24]. 
Further elaboration on this approach resulted in the first fully differential program to 
calculate an NNLO cross section [26, 28] and has lead to differential NNLO results for a 
number of processes meanwhile [29,31,38,41]. 



2 Basic concepts 

To introduce the subject, let us look at the simple example of a two-dimensional param- 
eter integral of the following form: 



dx 



dyx~ 1 ^y^(x+{l-x)y) \ (1) 



The integral contains a singular region where x and y vanish simultaneously, i.e. the 
singularities in x and y are overlapping. Our aim is to factorise the singularities for 
x — ► and y — > 0. Therefore we divide the integration range into two sectors where x 
and y are ordered (see Fig. [T]) 



l /.l 
dx 



dyx~^ ae y' be (x + (1 - x) y) 1 [@{x - y) + 0(y - x)] . 

(i) (2) 

Now we substitute y = x t in sector (1) and x — y t in sector (2) to remap the integration 
range to the unit square and obtain 

I = 



dxx- 1 '^ 6 dtr bf - (l + (1 - x) tj 



3 



y 



y 



Figure 1: Sector decomposition schematically. 

+ £ dyy- 1 -^* £ dtt- 1 - a *(l + (l-y)ty 1 , (2) 

We observe that the singularities are now factorised such that they can be read off 
from the powers of simple monomials in the integration variables, while the polynomial 
denominator goes to a constant if the integration variables approach zero. The same 
concept will be applied to TV-dimensional parameter integrals over polynomials raised 
to some power, where the procedure in general has to be iterated to achieve complete 
factorisation. 



3 The algorithm for multi-loop integrals 

3.1 Feynman parameter integrals 

A general Feynman graph Gf 1 '"^ 11 in D dimensions at L loops with N propagators 
and R loop momenta in the numerator, where the propagators can have arbitrary, not 
necessarily integer powers Vj, has the following representation in momentum space: 



G 



Hx—HR _ 



N 



1=1 



UPp({k}Ap},m]) 

3=1 



4-D 



— d u h , Pj{{k}, M, mj) = (qj - mj + iS) , 



(3) 



where the qj are linear combinations of external momenta pi and loop momenta k\. 
Introducing Feynman parameters according to 



r(jv„) 



nf=ir(^ 



N 



r oo N N 

/ n^-r^-E 



where N v = Vj , leads to 



3=1 



1 



P, 



A',, 



(4) 



■■■Mk 
1\...Ir 



fvi . . . rvj 

M * Ft. 



/ n^r^-E^)/ 



d^Ki . . . d^KL 



L L 

5^ kjM ij k j -2^2kj-Q j + J + i6 

*,i=i 3=1 



-Aft/ 



(5) 
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where M is a L x L matrix containing Feynman parameters, Q is an L-dimensional vector 
composed of external momenta and Feynman parameters, and J contains kinematic in- 
variants and Feynman parameters. The factors of k^ in the numerator can be generated 
from G(R = 0) by partial differentiation with with respect to Qf* , where I S {1, . . . , L} 
denotes the / th component of the vector Q, corresponding to the I th loop momentum. 
Therefore it is convenient to define the double indices Ti = (I, fii(l)) , I £ {1, . . . , L}, i £ 
{1, . . . , R} denoting the i th Lorentz index, belonging to the I th loop momentum. 

To perform the integration over the loop momenta ki , we perform the following shift 
in order to obtain a quadratic form for the term in square brackets in eq. ([5]): 



k[ = h-vi , vi = ^ M u 1 Qi 



(6) 



After momentum integration one obtains 



Ix—Ir 



(-1)*' 

[R/2] 



N 



N 



*(!-$>) 



jrN v -LD/2- 



1=1 



E (-^ mr ( N » LD i 2 ^ [^ rl ® -9) (m) 

m=0 

U N„-(L+l)D/2-R 

X 



2m) 



Ti,...,Tt 



(7) 



where 



^(x) = det(M) 
W(x*) = det(M) 



E 



Qj M7 1 Qi-J-iS 



= UM~ 



l = Uv 



(8) 



and [i?/2] denotes the nearest integer less or equal to R/2. The expression 

[(M^ 1 ® ^^( m ) ^ (-R-2m)jri,...,r B s t anc [ s f or ^ ne sum over a n different combinations of 

R double-indices distributed to m metric tensors and (i? — 2m) vectors Z. The above 
expression is well known [14,72-76], but an example to illustrate the distribution of 
indices may be helpful, so let us consider a two-loop integral where the fci-integral is a 
rank two tensor integral and the /^-integral is of rank one: 



^112 



d°K 1 d 1J K2 



JL/*1 

nj^ fv-^ nj*2 



N 



3=1 



11.7=1 j=i 1=1 



xi ) 



(9) 



{ 



7/JV„-3D/2-3 



j/N v -3D/2-3 

D) — rr — =— - x 



J" 



5 



Mf/ fit"" g" v 



The functions U and J 7 also can be constructed from the topology of the corresponding 
Feynman graph as follows [74,77,78]. Cutting L lines of a given connected L-loop graph 
such that it becomes a connected tree graph T defines a chord C(T) as being the set of 
lines not belonging to this tree. The Feynman parameters associated with each chord 
define a monomial of degree L. The set of all such trees (or 1-trees) is denoted by T\ . The 
1-trees T G T\ define hi as being the sum over all monomials corresponding to a chord 
C(T E T{). Cutting one more line of a 1-tree leads to two disconnected trees, or a 2-tree 
T. T2 is the set of all such 2-trees. The corresponding chords define monomials of degree 
L + 1. Each 2-tree of a graph corresponds to a cut defined by cutting the lines which 
connected the two now disconnected trees in the original graph. The momentum flow 
through the lines of such a cut defines a Lorentz invariant Sf = QZjecut(T ) Pj) 2 - The 
function !Fq is the sum over all such monomials times minus the corresponding invariant. 
For a diagram with massless propagators, T = J-q. If massive internal lines are present, 
T gets an additional term as follows: 

= sin*;]' 

TeTi j£C(T) 

= e [ n x ] ( 



Srp) , 



TeT 2 iec(T) 



N 

T{x) = Mx)+U{x)J2 x 3 m2 ■ (10) 

3=1 

U is a positive semi-definite function. Its vanishing is related to the UV subdivergences of 
the graph. Overall UV divergences, if present, will always be contained in the prefactor 
T(N U — LD/2). In the region where all invariants Sf are negative, which we will call 
the Euclidean region in the following, T is also a positive semi-definite function of the 
Feynman parameters Xj. Its vanishing does not necessarily lead to an IR singularity. 
Only if some of the invariants are zero, for example if some of the external momenta 
are light-like, the vanishing of T may induce an IR divergence. Thus it depends on the 
kinematics and not only on the topology (like in the UV case) whether a zero of T leads 
to a divergence or not. This fact makes it much harder to formulate general theorems 
for the subtraction of IR singularities of multi-loop Feynman graphs. The necessary (but 
not sufficient) conditions for an IR divergence are given by the Landau equations [79-81], 
which, in parameter space, simply mean that the necessary condition T = for an IR 
divergence can only be fulfilled if some of the parameters Xi go to zero, provided that all 
kinematic invariants Sf are negative. 

As can be seen from Eq. ([7]), the difference between scalar and tensor integrals is, 
once the Lorentz structure is extracted, given by the fact that there are polynomials 
of Feynman parameters in the numerator. These polynomials can simply be included 
into the sector decomposition procedure, thus treating tensor integrals directly without 
reduction to scalar integrals. 
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However, there is yet another possibility: Any tensor integral can be expressed in 
terms of scalar integrals in shifted dimensions, with some of the propagator powers 
different from unity, as has been shown in [75,78]. As our propagators can have arbitrary 
powers i/j, and the dimension D is a free parameter, this is a viable alternative. 

3.2 Iterated sector decomposition 

For less trivial examples than the one given in section [2| the singularities will not be 
factorised already after a single step of sector decomposition. An algorithm how to 
iterate this procedure is described below. 

Our starting point is a function of the form of Eq. ([7]) . As the basic algorithm is the 
same for tensor integrals, we will consider R = here for ease of notation, i.e. 

= ^ „v ,,„., ./ 1 ^ -X>> rw^m ■ 



3 = 1 1=1 



(11) 



Part I Generation of primary sectors 



We split the integration domain into N parts and eliminate the ^-distribution in such 
a way that the remaining integrations are from to 1. To this end we decompose the 
integration range into N sectors, where in each sector I, xi is largest (note that the 
remaining Xj^i are not further ordered): 



N „oo N 



pOO iV poo iv 

/ d N x = J2 d N xT[e{x l > Xj ) . (12) 
Jo ,_i Jo 



i=i " u j=i 



The (9-function is defined as 

0(x >y) = 



1 if x > y is true 
otherwise. 




The integral is now split into N domains corresponding to N integrals Gi from which 
we extract a common factor: G = (—l) Nv T(N l , — LD/2) 53;=i l n the integrals Gi 
we substitute 

for j < I 

for j = I (13) 
l for j > I 

and then integrate out xi using the ^-distribution. As U, T are homogeneous of degree 
L,L + 1, respectively, and xi factorises completely, we have U{x) — > liiit ) xf and T{x) —> 
Ti(t ) x[ +1 and thus, using J dxi/xi 6(1 — xi(l + ^2^=1 = 1) we obtain 

n .7=1 l J J 
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Note that the singular behaviour leading to 1/e-poles still comes from regions where a 
set of parameters {U} goes to zero. This feature would be lost if the (^-distribution was 
integrated out in a different way, since this would produce poles at upper limits of the 
parameter integral as well. The N generated sectors will be called primary sectors in 
the following. 

Part II Iteration 

Starting from Eq. (|14|) we repeat the following steps until a complete separation of 
overlapping regions is achieved. 

II. 1: Determine a minimal set of parameters, say S = {t ai , . . . ,t Qr }, such that Ui, 
respectively Ti, vanish if the parameters of S are set to zero. S is in general not 
unique, and there is no general prescription which defines what set to choose in 
order to achieve a minimal number of iterations. Strategies to choose S such that 
the algorithm is guaranteed to stop are given in [18] . Using these strategies however 
in general leads to a larger number of iterations than heuristic strategies to avoid 
infinite loops, described in more detail below. 

II. 2: Decompose the corresponding r-cube into r subsectors by decomposing unity ac- 
cording to 

r r r 

n *(i > *«, >o)=eii ^ *«* ^ °) • ( i5 ) 

j=l k=l i=i 

j^k 

II. 3: Remap the variables to the unit hypercube in each new subsector by the substitu- 
tion 

This gives a Jacobian factor of i^T 1 - By construction t ak factorises from at least 
one of the functions Ui , J~i . The resulting subsector integrals have the general form 

1 (N-l \ 7J N v -(L + l)D/2 

n *j i u- LD/ > > *=i,...,r. (it) 

For each subsector the above steps have to be repeated as long as a set S can be found 
such that one of the functions Ui... or Ti,„ vanishes if the elements of S are set to zero. 
This way new subsectors are created in each subsector of the previous iteration, resulting 
in a tree-like structure after a certain number of iterations. The iteration stops if the 
functions Uik 1 k 2 ... or ^ikik 2 ... contain a constant term, i.e. if they are of the form 

Ui klk2 ... = l + u{t) (18) 

Flkik 2 ... = so + ^2(-8p)fp(£) , 

P 

where u(t ) and fp(t ) are polynomials in the variables tj (without a constant term), and 
S/3 are kinematic invariants defined by the cuts of the diagram as explained above, or 
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(minus) internal masses. Thus, after a certain number of iterations, each integral G/ is 
split into a certain number, say a, of subsector integrals. We can replace the multi-index 
fci&2 • • ■ stemming from the subsector decomposition by a single index which just counts 
the number of generated subsectors. The subsector integrals are exactly of the same 
form as in Eq. (J7J), with the difference that the index k now runs from 1 to a, the total 
number of produced subsectors in each primary sector. 

Evidently the singular behaviour of the integrand now can be read off directly from 
the exponents aj , bj for a given subsector integral. As the singular behaviour is manifestly 
non-overlapping now, it is straightforward to define subtractions. 



Part III Extraction of the poles 

The subtraction of the poles can be done implicitly by expanding the singular factors 
into distributions, or explicitly by direct integration over the singular factors. In any 
case, the following procedure has to be worked through for each variable ij=i,...,jv— 1 and 
each subsector integrand: 

• Let us consider Eq. (|17p for a particular tj, i.e. let us focus on 

1 

J, = / dtjtp-Wlfaitv&e), (19) 



where T = U^"~^ L+1 /J 7 ^." LD ^ 2 in a particular subsector. If aj > — 1, the 
integration does not lead to an e-pole. In this case no subtraction is needed and 
one can go to the next variable tj+i. If aj < — 1, one expands 1{tj, {U^j}, e) into 
a Taylor series around tj =0: 

Kl-i p 



&H0,{t*iU) - d?I(t J ,{t l ^},e)/dt P 1 . (20) 



f,=0 



• Now the pole part can be extracted easily, and one obtains 



'> - ^ ^U-^^'^ h^^ ■ (21) 

p— o 

By construction, the integral containing the remainder term R(t, e) does not pro- 
duce poles in e upon ^-integration anymore. For aj — —1, which is the generic 
case for renormalisable theories (logarithmic divergence), this simply amounts to 



h 



i 

^(0, {t t ^}, e)+J dtj tj 1 -^ (ifaiUtf}, e) - 1,(0, {U^}, ej) , 



which is equivalent to applying the "plus prescription" [82] (see eq. (|45)l ). except 
that the integrations over the singular factors have been carried out explicitly. 
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Since, as long as j < N — 1, the expression (f2"Tj) still contains an overall factor 
ij+i e 3+1 , it is of the same form as (fT£))) for j — > j + 1 and the same steps as 
above can be applied. 

After iV — 1 steps all poles are extracted, such that the resulting expression can be 
expanded in e. This defines a Laurent series in e with coefficients Ca-,m for each of the 
a{l) subsector integrals Gik- Since each loop can contribute at most one soft and collinear 
1/e 2 term, the highest possible infrared pole of an loop graph is l/e 2L . Expanding to 
order e r , one has 

2L „ N a(l) 

G ^= E ^ + 0(e r+1 ), G={-l) N »T(N v -LD/2)Y,Y, G u- ^ 

m= — r 1=1 fc=l 

Following the steps outlined above one has generated a regular integral representation of 
the coefficients C^-.m, consisting of (N — 1 — to) -dimensional finite integrals over param- 
eters tj . We recall that T was non-negative in the Euclidean region where all invariants 
are negative (see eqs. (|10I18[1 ). such that the numerical integrations over the finite pa- 
rameter integrals are straightforward in this region. In principle, it is also possible to do 
at least part of these parameter integrals analytically, but in most applications such an 
analytical approach reaches its limits very quickly. 



Avoiding infinite recursion 

As mentioned already, the choice of the set S = {t ai , . . . ,t ar } which makes U respec- 
tively T vanish for t a — > is in general not unique. The structure of the function hi (see 
cq. (|10p ) is such that its decomposition will always terminate after L iterations for an 
L-loop integral. For the function T ', the structure depends on the masses and kinematic 
invariants involved. Although one could follow one of the mathematical strategies given 
in [18] to ensure the iteration terminates, this is not the most efficient method for practi- 
cal purposes, as these strategies typically generate a large number of subsectors. Another 
possibility, adopted in [25], is to choose the set S randomly, such that eventually a set 
will be selected which does not lead to infinite recursion. However, it is more efficient to 
use some heuristic rules which, in all applications to multi-loop diagrams considered so 
far by the author, lead to a terminating decomposition procedure. 

Let us first illustrate the problem by a simple example: Consider the function 

f{xx,X2,x 3 ) = x\ +x%x 5 , (23) 

and suppose we choose S = {1,3}. The replacement x% — x^tx in the subsector as- 
sociated with 9(x 3 — xi) leads to / = X3 (13 tf + x\). Choosing now S = {2,3} and 
substituting x-i = x^t 2 in the corresponding subsector brings us back to the original 
functional form, so we generate an infinite recursion for the above choices of S. In this 
simple example we can see immediately that the choice S = {1,2} does not lead to this 
problem. 

For multi-loop integrals, we can use the following facts as a guideline to choose con- 
venient sets S: We first note that an infinite recursion does not occur for functions which 
are linear in each variable. The function T, before the iterated decompositions, is a poly- 
nomial of maximal degree two in each individual Feynman parameter, where quadratic 
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N 

parameters only occur if massive propagators are present, due to the term U (x) ^ Xjirij 

contained in T (see eq. (|10|) '). Therefore, a simple extra rule for diagrams with internal 
masses can be added to the procedure: Before entering the iteration, determine the set 
Sm of labels belonging to massive propagators and use this set for a first sector decom- 
position (even if it does not lead to T — upon setting the elements of Sm to zero). 
This produces a form where in each subsector, one of the quadratic powers is reduced 
by one, such that self-similarity to the original form cannot be generated anymore. In 
the course of the iterations, quadratic or higher powers will be generated unavoidably, 
such that a form which may lead to infinite recursion can occur at some point. In this 
case it has proven useful to choose, if existent, a set S containing the maximal number of 
variables occurring with the same power. Certainly, these are only heuristic rules, which 
however worked well in a multitude of practical applications. 

3.3 Examples 

3.3.1 Planar double box with one off-shell leg 

Two-loop box diagrams with one off-shell leg are master integrals entering for example 
the two-loop QCD matrix elements for e + e~ — ► 3 jets at NNLO [47]. Numerical results 
were first given in [16] and served as important benchmarks for the analytical calculations 
of Refs. [45,46]. 
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Figure 2: The planar double-box with leg 4 off-shell. 

For the planar double box with p\ = p\ = p\ = 0,p| ^ shown in Fig. [3J the 
functions IA and T are given by 

U = Xi 23 X 567 + X 4 Xi23567 

T = (-si 2 )(x 2 x 3 a;45 6 7 + £5X6X1234 + x 2 x 4 x 6 + x 3 x 4 x 5 ) 

+ (-S23)xiX 4 X7 + (-P4)X7(X 2 X 4 + X5X1234) , (24) 

where x ijk ... = x t + Xj + x k + . . . and s tj = {p t + pj) 2 . 

Iterated sector decomposition produces 197 sectors. As the off-shell leg regulates 
some of the singularities which would be present in the planar double box with all legs 
on-shell, the number of produced subsectors is lower than for the on-shell planar double 
box (282 subsectors) . The result for two Euclidean points is given in Table [U where an 
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overall factor of T(l + e) 2 has been extracted and the integral is defined atQ 



d D h fd D k 2 

DB mi = / — =- / —s- X 



k\ k\ (h + Pl y (h + Pl + P2 y (fc 2 + Pl + P2 y (k 2 - Pi y (h - k 2 y 
r(i + *) a (5 + 5 + 5 + ? + ^ ■ (25) 



The computing time for the given precision, which is better than 0.3% for the finite 

~ls^rs 23 ,s 13 ,pi) (-1/3,-1/3,-1/3,-1) (-1/2, -1/3, -1/6, ^lf 

Pi -26.9997±0.00049 -11.9998±0.0002 

P 3 -118. 651±0. 0037 -43.0010±0.0027 

P 2 -239.646±0.0347 -58.6686±0.0160 

Pi -305.823±0.1835 -20.7692±0.0560 

Pa -162.537±0.435 +98.191±0.289 



Table 1: Numerical results for the pole coefficients of the planar double-box with one leg 
off-shell. An overall prefactor of T 2 (l + e) has been extracted. 

part and better than 0.1% for the pole coefficients, was about 2hrs on 3.0 GHz Intel 
Xeon processors. To obtain a precision of only 1% in the finite part takes about 30 
minutes. The numerical evaluation has been done for each primary sector separately 
and the errors have been added in quadrature. The independent treatment of each 
primary sector allows to split the problem into smaller subparts which can be evaluated 
in parallel, such that the overall computing time is determined by the primary sector 
with the most complicated singularity structure. Further, symmetries of the diagram can 
serve as a check, as the results for the corresponding primary sectors should be identical. 
On the other hand, if large cancellations between different primary sectors are observed, 
summing over the primary sectors before the numerical integration is the better option. 



3.3.2 Three-loop vertex diagram 

As a more complicated example, let us consider the diagram shown in Fig. [3l entering 
the calculation of massless three-loop form factors [55]. It is given by 

J {2rrf J {2tt) D J {2rrf 

1 

(fcl + pi) 2 k\ k\ (k 2 + Pl ) 2 (fc! + fc 3 ) 2 Wl + h + qf (k 2 - fci) 2 (k 2 + fc 3 ) 2 ' 

where q = p\ + p 2 is the incoming momentum, and again, an infinitesimal imaginary 
part +i8 in the propagators is understood. Iterated sector decomposition produced 684 
sectors. The result is given in Table [5J where an overall prefactor of i Sp (— q 2 — i5)~ 2 ~ 3e 

4 An infinitesimal imaginary part +iS in the propagators is understood, and we use jj, = 1. 



12 



Figure 3: A non-planar 3- loop vertex diagram, called Ag. 

with 1/Sr = (47r) D / 2 r(l - e), has been extracted: 

A s = iS$ (-q 2 - l 5)-^ + ^ + P + eP, E ) , (27) 

and the Pj are the coefficients given in Table [3J The computing time up to order e was 
about 4hrs on 3.0 GHz processors. Computing only the pole coefficients and the finite 
part took about 1 hour. 





numerical 


analytic 




3.20553 ±0.00011 


3.2054850751 


Pi 


8.42310 ±0.00146 


8.4222653365 


Po 


27.885 ± 0.039 


27.852843117 


Pe 


-50.246 ±0.129 


-50.283167385 



Table 2: Numerical results for the Laurent-expansion of the 3- loop vertex diagram shown 
in Fig. [3l The analytic result can be found in [55]. 

4 Sector decomposition for infrared divergent real ra- 
diation integrals 

In order to calculate cross sections at higher orders in perturbation theory, there are 
in general not only virtual corrections, but also corrections from real radiation to be 
taken into account. At next-to-leading order, we only have two types of contributions: 
the purely virtual (one-loop) corrections, and the real radiation of one additional par- 
ticle, which may be either theoretically or experimentally unresolved. "Theoretically 
unresolved" denotes the collinear branching of massless particles or the emission of soft 
photons or gluons, which leads to infrared singularities appearing as poles in 1/e in 
dimensional regularisation when integrated over the D-dimensional phase space. Exper- 
imentally unresolved particles do not lead to a 1/e-singularity. They are defined by a 
so-called "measurement function" defining the physical observable, which in most cases is 
a subroutine in a numerical program rather than an analytic function. For example, two 
particles which are clustered into a single jet by a certain jet algorithm are considered 
as experimentally unresolved. 

At NNLO, one generally has to deal with three building blocks making up the full 
cross section: two-loop (and one-loop squared) virtual corrections, one-loop virtual cor- 
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rections combined with single unresolved real radiation, and doubly unresolved real ra- 
diation. 

As any _D-dimensional phase space integral can be transformed to a dimensionally 
regulated multi-parameter integral over the unit hypercube, the singularities stemming 
from the (theoretically) unresolved real radiation are amenable to sector decomposition 
applied to phase space integrals over the corresponding squared matrix elements. What 
matters here are the denominators of the matrix elements for different processes, which 
have a generic form, and therefore allow for the setup of a general framework. 



4.1 Phase space integrals in D dimensions 

The phase space integral in D dimensions for a generic process Q — + p\ + . . . + pn can 
be written as 

. ~ N N 

J J j=l i=l 

where 5 + (p 2 — to 2 ) = S(p 2 — m 2 )0(p'°)). Using 



d u p 3 5+{pj-mj) = — j d»~% 



E 3 = x/P? +m2 j 

for j = 1, . . . , N — 1 and eliminating p^ by momentum conservation, one obtains 



d^, = 



N-DIN-1) fN- 1 Q/nl . N ~! 



■n - / 11" fi E . u Z^M) - ""N 

J .7=1 3 i=l 



To proceed, one has to choose a certain parametrisation for the phase space integration 
variables and work out the integration limits confining the integration range to the 
physical region. The scattering case Q = p a + Pb ^ N — 1 particles differs from a decay 
1 — > N particles by the fact that in the center-of-mass frame of the incident particles, it 
contains a preferred direction given by the beam axis p a = —pt- Finding the appropriate 
phase space integration variables which are optimally adapted to the kinematic situation 
at hand can simplify the calculation considerably. This is even more true if sector 
decomposition is used to isolate the infrared singularities: a convenient parametrisation 
will be one where the maximal number of potentially singular denominators of the matrix 
element naturally factorises, thus limiting the number of terms produced by iterated 
decompositions. In fact, it turns out to be useful to divide the matrix element into 
different "topologies" , according to their denominator structure, and use several phase 
space parametrisations, each being optimal for a certain class of topologies. 

A multi-particle phase space is most conveniently described as a convolution of phase 
spaces of lower multiplicity. For example, a process like the one in Fig. [4] suggests a 
phase space parametrisation which is a convolution of a phase space for a 1 — > 4 decay 
followed by a 1 — > 3 decay and a 1 — > 2 splitting. For a process involving soft radiation 
off massive fermions, it is convenient to choose a parametrisation where the energy of the 
particle which can become soft is an integration variable. Useful examples of different 
parametrisations can be found e.g. in [28,83,84]. Here, in order to illustrate some generic 
features of the method, we will first derive a phase space parametrisation in terms of 
double invariants Sij = (pi +pj) 2 . 
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Figure 4: Example of a cascade decay corresponding to an iterative construction of a 
multi-particle phase space. 



As a pedagogical example we will consider the massless case, pj = 0. Let use choose 
a 1 — > 4 process and consider a special frame where 

Q = (E,Q( D -V) 

pi = Ex{l^ D - 2 \l) 

P2 = £2(l,O (c - 3) ,sin0i,cos0i) 

p 3 = E 3 (l,0 (D - 4 \sin6bsin6> 3 ,sin(9 2 cos6> 3 ,cos6>2) 

Pa = Q - Pi - Pi - P3 , (29) 



which leads to 



1 



<i$i^4 = -{2tt) a ~ zd dE 1 dE 2 dE 3 d9 1 d9 2 d9 3 [E 1 E 2 E 3 sin9 1 sin9 2 } D - 3 sm9 3 ~ 4: 
8 

dn D _ 2 dn D _ 3 dn D ^<d{E 1 )Q{E 2 ) e(E 3 )e(E ~e 1 -e 2 - e 3 ) 

5{E 2 ~ 2E(E 1 +E 2 + E 3 ) + 2{ Pl -p 2 + Pl -p 3 +p 2 - p 3 )) , (30) 
where 

dClr,-! = / d6»i / d6 2 sm6 2 ... d9 D - 1 (sm9 D - 1 ) D -' 2 - 



r(#) 



Now we map the angle and energy variables to the double invariants Sy as integration 
variables, using the Jacobian 

det(J) = det [ f^ S ~l ) = 64 E 3 E 2 E 2 E 2 sin6> 2 sinfl 2 sin6> 3 . (31) 



d{Ei,6 3 ) 

The Jacobian in combination with terms already present in (|30|) can be written as the 
determinant A4 of the Gram matrix = 2pi -pj. This determinant can be expressed 
by the Kallen function X(x,y, z) = x 2 + y 2 + z 2 — 2xy — 2yz — 2xz as 

A 4 = A(si2 s 34 , S13 s 2 4, S14 S23) = - (4EEi E 2 E 3 sin#i sinfo sin^) 2 . (32) 

We see that A4 has to be negative semi-definite. With the dimensionless variables 

yi = s 12 /Q 2 , y 2 = s 13 /Q 2 , y 3 = s 23 /Q 2 , y 4 = s u /Q 2 , ys = s 24 /Q 2 , 2/6 = s 34 /Q 2 

(33) 
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and A 4 = A(yiy 6 ,y 2 2/5,?/32/4) we obtain finally 



(2n)*- 3D (Q 2 ) 3D/2 - 4 2- 2D+1 dSl D _ 2 dQ D . 3 dtl D _ A 



3=1 



9(-A 4 ) [-k^ D -^H{l - J2 Vi) ■ (34) 

3=1 



For N > 5 we have to distinguish if we are in D dimensions or in four dimensions. 
In D dimensions, the same procedure as above can in principle be applied. The four- 
dimensional case is complicated by the fact that the Gram determinant A^r vanishes for 
N > 4. In this case the phase space can be expressed in terms of the Kallen function of 
invariants built from four independent momenta and additional constraints [83], but in 
practice it is more useful to build it up iteratively as described above. 



4.2 Special features of sector decomposition for real radiation 

We see that expression (|34l) has a high symmetry in the invariants yj. To proceed in 
a way analogous to the treatment of loop integrals, we could now do a "primary sector 
decomposition" to integrate out the <5-function as explained in section 13.21 This would 
lead to n primary sectors, where n is the number of two-particle invariants s^, i.e. n = 6 
in the example above. All invariants are treated on equal footing in this step. The 
primary sector decomposition is very useful in the case of loop integrals, mainly for the 
following reason: it preserves the feature that singularities only occur at special points 
at the boundary of parameter space: they occur only if y^ , . . . , y, w — for a subset 
{it . . .i r } of {1 . • .n}. In other words, in the case of loop integrals in the Euclidean 
region, no singularities can occur for yi — > 1 or in the interior of parameter space, and by 
primary sector decomposition the ^-constraint is integrated out without destroying this 
feature. In addition, the integration limits from zero to one for all remaining variables 
are guaranteed without further transformations. 

In the case of real radiation, the situation is different, because we are forced to stay 
within the physically allowed region. In the parametrisation above, this is reflected 
by the fact that after integrating out the constraint 5(1 — X^i^) fr° m momentum 
conservation, we still have the constraint 0(— A4). Solving the equation A4 = for say, 
yk, we obtain the solution y k = (y/x ± \/~z) 2 , so y^ = whenever x — z. For example, 
for k — 3, we have 

yf = (VV2V5 ± ^Viyaf/Vi ■ (35) 

The substitution y 3 = y^ + t 3 (yf — y$ ) in order to remap all integrations to the unit 
interval will lead to a complicated structure of those denominators in the matrix element 
which contain j/3. In fact, we see that l/j/3 will develop a singularity if £3 = and y^ 
simultaneously, i.e. whenever £3 = and 2/22/5 = 2/12/6- Thus we found two properties 
which did not occur in the case of loop integrals: 

1. Square-root terms appear naturally when solving the phase space constraints. Such 
terms are potentially dangerous as they may destroy the polynomial structure 
which is a prerequisite for the sector decomposition, leading to expressions like 
g(x, y) = a + y — y/a 2 + x 2 , where a is a constant. However, it is obvious that such 
terms can be easily transformed into a form with the required behaviour under 
rescaling of the variables. 
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2. After having mapped the phase space integration limits to the unit hypercube, 
singularities can occur for a manifold which (partly) lies inside the phase space 
integration region. 

How these singularities can be remapped to the boundaries will be shown in a specific 
example below and discussed more generally in Section [4~4l Here we would like to point 
out that we cannot solve the constraint 0(— A4) for the same variable yk in each primary 
sector, because in primary sector k, yk has been eliminated. Therefore a judicious choice 
of yk — to be an invariant which occurs only in very few or no denominators of the 
complete matrix element — would still lead to complicated denominators in primary 
sector k, where the constraint had to be solved for yi-tk- For this reason it is advisable 
not to use primary sector decomposition in the case of complicated matrix elements for 
real radiation. 





loop integrals 


phase space integrals 


parametrisation 


hi functions in terms of 


many different options, should 




Feynman parameters fixed 


be adapted to topologies 


primary sector dec. 


very convenient 


not recommended 


singularity structure 


in Euclidean region only 


singularities inside integration 




endpoint singularities 


region generic 



Table 3: Main differences to loop integrals in the sector decomposition procedure for 
phase space integrals over IR divergent real radiation 

To choose a parametrisation which is adapted to the denominator structure of the 
problem, one can follow the idea of iterative splittings outlined above. Matrix elements 
involving massless particles contain invariants of the form Si 1 ...i n — {pn + ■ ■ • +Pi„) 2 with 
11 > 2 in the denominator. For example, the four-particle cut of the diagram in Fig. [5] 




-wwvwu N vwvwv^ 



Figure 5: Example of a four-particle cut. 
contains an integral of the form [23] 



J 4 = - /nd yi e(-A 4 )(-A 4 )- 1/2 -^(i-& 

^~ J 1 A. 1 



(2/1 + 2/5) (2/2 + ye) - 2/32/4 



,_, 2/2(2/2 +2/4 + 2/e) 2 

J ~ 

(36) 

In this case, it is suggestive to introduce the triple invariant S134/Q 2 = 2/2+2/4 + 2/6 as 
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a genuine phase space variable, such that this denominator factorises immediately. This 
example will be worked out in detail in the following section. 

Obviously, it is advantageous to use triple invariants as phase space integration vari- 
ables if the amplitude contains a splitting of one particle into three final state particles, 
double invariants if the amplitude contains several 1 — > 2 splittings, etc. Therefore the 
choice of parametrisation is most conveniently done on a topology basis, i.e. different 
parametrisations are applied to certain classes of denominator structures, as already 
mentioned above. As the full matrix element contains interferences of amplitudes of 
different "splitting history" , it is in general impossible to achieve a factorised form for 
all denominators. However, minimising the number of decompositions by convenient 
parametrisations is vital to limit the size of the expressions produced by iterated sector 
decomposition. 

The main differences to loop integrals in the sector decomposition procedure for phase 
space integrals are again summarised in Table [3] 

4.3 Example of a four-particle final state 

To explain the concept, we go back to the example of the previous section, the massless 
1 — > 4 phase space, and topologies containing S134 or S234 in the denominator. In order 
to achieve a convenient parametrisation, we first multiply eq. ([51)1 by 

1 = J dx A S(x 4 - j/2 - 2/4 - 2/6 ) J dx 5 S(x 5 - y 3 - y 5 - y 6 ) 

and eliminate y±, j/4, j/g using the S- functions (see eq. Q33p for the definition of the scaled 
invariants). Then we solve the constraint (— A4) > for j/3 = S23/Q 2 and substitute 
j/3 — > j/g" + i 3 (2/3" — 2/3"). For the remaining variables we substitute 

X4 = ti 

x 5 = h + 1 6 (l - fi) (l - U) 
yi = h(l-t 2 )t4 

y G = hU (37) 
to arrive at the following form for the phase space 

^$1^4 = (2n) i - 3D {Q 2 ) 3D ^- 4 2- 2D+1 dn D ^ 2 dn D _ 3 dn D _ 4 (38) 

. . . dt 5 [(i - i x )t 4 (i - h)] D - 3 [tMi - hW - h)}^ 1 [«3(i - t 3 )] (D - 5)/2 . 

The expression for J4 in eq. (|36[) then becomes 

dh... dt 5 [h(l - t 2 )]- 1 -% 1 ~^(l - h) 1 - 2 ^! - U) 2 ' 2 ' 

[*2* B (i - h)r e Mi ~ h)]- 1/2 - e {(h(i - 1 2 ) + 1 2 ) - fe(t)) (39) 



1 

4 

7T 




y 3 (t) = hh + h(l - ta)(l - t 5 ) - 2(1 - 2*3)^1*2*5(1 - *2)(1 - h) 

= lft(t)/(l-t4). 
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We see that in this paramctrisation, the denominators in J4 are factorising completely. 
However, other denominators in the full matrix element will in general contain §3(1) in 
the denominator. In this case it is convenient to shuffle the square-root terms to the 
numerator by the following non- linear transformation [22]: We substitute 

2/3 + z 3 (2/3 - 2/3 ) VS + z 3 (2/3 - % ) 

The Jacobian therefore cancels one factor of 2/3 in the denominator: 

d*3 _ 2/3" 2/3" 2/3 



d ^3 [y 3 + z 3 (y+ - y 3 )] 2 % + z 3 (y+ - y 3 ) 

leading to 

/ dt 3 [t 3 (l-t 3 )}^ = [ dt 3 [i 3 (l - t 3 )] 

Jo 2/3 (*) Jo 



(41) 



1 



% + *3 (2/3 - 2/3 ) 

= [dz 3 [z 3 (l-z 3 )y+y 3 ] £ ^[y 3 +z 3 (y+-y^)} i - D . (42) 
Jo 

This way the square-roots in the denominator are eliminated and the limits t 3 — ► and 
7/3" — > are decoupled, but note that instead of l/y 3 (t) we now have a factor 

[2/ 3 + 2/3l ^ = [(1 - U) (Ml - *2)(1 - h) - t 2 t 5 )] D - 5 = (1 - t 4 ) D - 5 [f(t u t 2 , t 5 )f- 5 . 

The factor (1 — t^) ^ 5 will be combined with phase space factors and is of endpoint- 
type anyway, but there are singularities which now occur on a manifold defined by 
f(ti,t 2 ,t 5 ) = 0. In the case at hand they are easily remapped to the boundaries by 
splitting e.g. the ^-integration region at 

= *i(l-*5) (43) 

*5 + *l(l-t5) 



and substituting 



t 2 = t% u 2 for t 2 < t 2 , 
t 2 = 1 - (1 - t° 2 ) u 2 for t 2 > t° 2 (44) 



to obtain again integrals from zero to one. 



4.4 Possible types of singularities and their treatment 

As we have seen in the previous section, we have to deal with the following types of 
singularities: 

• endpoint singularities 

• singularities on a manifold not confined to the boundaries of phase space, more 
precisely the boundaries after having solved all constraints and remapped the in- 
tegrations to the unit hypercube. 
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Endpoint singularities, if not factorising from the start, are easily extracted by the sector 
decomposition algorithm. It should be mentioned however that, if we do not use primary 
sector decomposition, endpoint singularities can occur not only if an integration variable 
goes to zero, but also at the upper integration boundary (which is equal to one, after 
appropriate remapping) . In order to apply the algorithm described in section 13.21 we 
should remap the singularities for tk — > 1 such that they occur at the origin only. As 
some variables can cause singularities at zero and one, a transformation tk — > 1 — tk is 
not recommended. Instead, we split the integration range at 1/2: After the split 



and the substitution tk = Uk/2 in (a) and tk = 1 — Uk/2 in (&), all endpoint singularities 
occur at Uk — > only. The disadvantage of such splittings is the fact that we end up 
with 2™ integrals after n splittings, but in practice, considering the physically possible 
singular limits, some of the integration variables clearly never will lead to a singularity 
when approaching one, and therefore do not require such a splitting. 

Concerning the singularities at the interior of phase space, the recipe is less simple. 
However, for A r <5inal->iVora2^A r -l phase space, it is easy to see that 
they can always be remapped to the phase space boundaries. Quite in general, the 
boundaries of the physical region in the space of invariants follow from the momentum 
conserving (^-function and the Gram determinants Ajv = A(pi, . . . ,pn) — det(p$ • pj). 
As the Gram matrices are symmetric, the determinants will be polynomials of maximal 
degree two in each invariant Sij. Masses do not alter this argument, and one can show 
that always A3 > and A4 < [84]. The case N = 3 is trivial, so let us first consider 
the case N = 4. Solving the constraint 8(— A4) for one of the invariants yk leads to 

= {\f a ~k± Vbk) 2 /ck, where the structure of a/-, bk, Ck is fixed by the fact that A4 < is 
a Kallen function: these terms must be linear in each invariant (see e.g. eq. (|35p ). After 
having performed substitutions of the type (f3"T)) to eliminate the momentum conserving 6- 
function, the linearity is not manifestly preserved, but as the singularity structure cannot 
change by these substitutions, some of the ti must always factorise, which guarantees 
that the condition A4 < imposing the phase space boundaries in the new variables can 
be solved for, say, the variable t® in such a way that i° is the ratio of two polynomials 
in the remaining parameters (see e.g. eq. (|43p ). therefore leading to structures amenable 
to sector decomposition. 

For A > 5, in a 1 - > iV or a 2 - *■ A — 1 phase space, additional constraints are 
present due the fact that, for 4-dimensional momenta, A^r = for N > 5. However, if 
the phase space is expressed in terms of a convolution of processes of lower multiplicity 
as explained above, the same reasoning for the remapping of singularities as in the N < 5 
case can be applied. 

Different phase space parametrisations are related by Lorentz transformations, there- 
fore it is sufficient to show this property for a particular parametrisation. This does not 
mean that all parametrisations actually do have the desired properties, it only states 
that a better parametrisation must exist where the remapping to a form more suitable 
for sector decomposition is possible. 
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One last point concerning different types of singularities should be made: In renor- 
malisable theories, the "physical" singularities are not worse than logarithmic, which 
means that the parameter integrals after sector decomposition should be of the form 
L dx x a+be f (x) , where a > — 1. If only the denominators are considered in a complex 
matrix element, terms with a < — 1 will occur. This type of spurious singularity will 
finally cancel with terms from the numerator, but the cancellation is not manifest if we 
leave the numerator symbolic throughout the whole procedure. Therefore it is advisable 
to include the numerator at the level of the e-expansion, at least for the parts where 
a < -1. 



4.5 Construction of a differential Monte Carlo program 

The isolation of infrared poles by sector decomposition is an algebraic procedure, leading 
to a set of finite functions for each pole coefficient as well as for the finite part. The finite 
functions have the form of parameter integrals over the unit interval and are therefore well 
suited for integration by Monte Carlo methods. If a full cross section beyond the leading 
order, composed of both real and virtual corrections, is to be calculated, the combination 
of the sector decomposition approach for the real radiation part with analytic results (if 
available) for the part involving loop corrections is certainly possible. In the case of 
NNLO corrrections, it is also advisable, as the fully numerical evaluation of two-loop 
integrals in combination with the phase space integration is in general rather slow, if 
viable at all. However, at NLO, examples of a fully numerical evaluation of complete 
processes based on the combination of sector decomposition with contour deformation 
do exist [58,59]. 

A hybrid approach can consist for example in the reduction of the phase space in- 
tegrals to cut master integrals, evaluating only the master integrals by sector decom- 
position [21]. Concerning the mixed one- loop times single unresolved real radiation 
part of NNLO calculations, its treatment so far always involved a reduction to master 
integrals [24,28,29], except in the very recent calculation of the 0(a 2 s ) corrections to 
semileptonic decay b — > clvi [41]. For these mixed real- virtual contributions, it can fur- 
ther be useful to do parts of the Feynman parameter integrations for the master integrals 
analytically, to obtain Hypergeometric functions where transformation formulas for the 
arguments [85,86] can be used if necessary to arrive at a more convenient form. Then 
one can proceed by applying sector decomposition to the integral representation of the 
Hypergeometric functions in combination with the phase space variables [29,87]. 

To obtain differential results, the combination of the output of the sector decompo- 
sition procedure with any infrared safe measurement function is possible, as has been 
first noted in [22]. The flexibility to do so is achieved by expanding the singular factors 
produced by the decomposition into plus-distributions, using the identity 



1 oo 

= _L six) + 



(«e)" 



IV. 



n=0 

where 

r- 1 /-i 



In" (a) 



dxf(x) [g(x)/x} + = I dx f(x) /(0) g(x). (45) 

x 







This is basically equivalent to the e-expanded form of eq. (f2"Tj) . the only difference being 
that, instead of integrating out the singular factors explicitly, the integrands are kept in 
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the form of distributions: dx x 1+K£ /(0, y) is written as J Q da; ^_ /(a;, instead 
of /(0, y). This allows for the combination with arbitrary functions f(x, y). 

The following features which are special to the sector decomposition approach as com- 
pared to Monte Carlo programs based on analytic subtraction terms should be pointed 
out: 

• The pole coefficients are only calculated numerically, such that the cancellation 
of poles between real, real-virtual (existing beyond NLO only) and purely virtual 
contributions can be verified only numerically. However, this is in general not a 
problem because the pole coefficients contain less integration variables and there- 
fore a high numerical precision can be achieved more easily than for the finite 
part. 

• The expansion into plus distributions cannot be done in complete isolation from 
the measurement function, because it has to be assured that the subtraction terms 
only come to action in phase space regions which are allowed by the measurement 
function. To illustrate this point, consider the simple one-dimensional example 
where the measurement function is just a step function Q(x — a), a > 0, and the 
"matrix element" after sector decomposition is given by f(x). If we expand into 
plus distributions and afterwards just multiply with our measurement function, we 
obtain 

/ dx f{x) - /(0) 6(x - a) = f (0) In a + / dx M . (46) 
Jo x J a x 

Clearly, the /(0) term stems from the subtraction of a singularity at x = 0, which 
is now killed by our measurement function anyway, such that inclusion of the /(0) 
term would lead to a wrong result. The correct way is of course to include the 
measurement function into the expression the plus distribution acts on. However, 
this does not mean that the e-expansions and subtractions have to be redone 
each time the measurement function is changed. It can be achieved by including 
symbolic functions in the e-expansion which are written to the numerical code with 
zero arguments (respectively the appropriate singular limit in the general case) if 
they correspond to subtraction terms. The symbolic functions can be specified 
later in a subroutine of the numerical program. 

• The subroutines defining jets, observables etc. will be based on the four-momenta 
of the particles involved in the scattering process. The four-momenta can easily 
be constructed from the original phase space integration variables. Before the 
decomposition, the phase space integration variables, let us call them s^ , have 
a certain functional form, Sij — Sy(ii, . . . , t n ). Performing now iterated sector 
decomposition will remap the parameters ti, in a different way in each subsec- 
tor of the decomposition tree, such that after iterated sector decomposition, the 
functional dependence of the original variables on the Monte Carlo integration pa- 
rameters t\, . . . , t n is different for each subsector. Of course it is easy to keep track 
of the remappings done in each sector, but the Monte Carlo program will consist 
of a sum of contributions from each subsector k, each one defining the functional 
form s\j\ti, ■ ■ ■ ,t n ) in a different way. This is not a problem in principle, but the 
complexity of NNLO matrix elements is already enormous, so multiplying the eval- 
uation time by the number of subsectors, which is of the order of several hundreds 
for an NNLO process, can lead to unacceptable CPU times. 
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As the subtractions done after sector decomposition are of the form 

m /(o) 



/ 

Jo 



dx ■ 



in each variable, which means that poles in each variable are locally subtracted, the 
method in general leads to expressions which have a good numerical behaviour. In 
fact, as even integrable singularities of the type f Q dxdy are decomposed, the 
expressions produced by iterated sector decomposition are of a form which is very 
convenient for numerical integrations. However, if the matrix elements to evalu- 
ate exceed a certain degree of complexity, there is a turnover where the advantage 
gained from the form of the individual functions is destroyed by the sheer number of 
functions to evaluate. This has been found for example in the attempt to calculate 
the full real corrections for e+e~ — > 3jets at NNLO using only sector decomposi- 
tion. The calculation of this process has recently been accomplished [88-91] using 
analytic "antenna" subtraction [92]. Due to the large number of massless parti- 
cles, the infrared structure is extremely complicated, and the number of antenna 
subtraction terms needed for the analytic subtraction of the poles is already quite 
large. Using sector decomposition leads to an unacceptable number of terms in this 
case. On the other hand, if massive particles are involved, the situation is com- 
pletely different: while analytic integrations of subtraction terms become nearly 
impossible for NNLO calculations with several mass scales, the infrared singularity 
structure is less complex in the presence of masses, such that the number of terms 
produced by sector decomposition will be moderate, and the mass dependence of 
the finite terms produced by the decomposition does not pose a problem for the 
numerical integration. 



5 Conclusions and outlook 

The method of sector decomposition is interesting from a more formal field theoretical 
point of view as well as for phenomenological applications. Within the context of dimen- 
sional regularisation, it offers a constructive scheme for the factorisation and subtraction 
of infrared poles to (in principle) all orders in perturbation theory, not only for individual 
integrals, but also for entire squared matrix elements. 

Quite in general, the method consists of two parts: the first is an algebraic one, where 
the singularities are isolated in terms of a Laurent series in e, the coefficients being finite 
parameter integrals. The second part consists in the evaluation of these parameter 
integrals, which in general is not possible analytically, and is therefore done numerically 
by Monte Carlo integration. Obviously the precision which can be achieved this way is 
intrinsically limited, compared to the evaluation of functions where all integrations have 
been performed analytically, or where deterministic numerical integration methods can 
be applied. However, in most practical applications considered so far, sufficient precision 
could be reached within a reasonable amount of integration time. 

Applications of the method to multi-loop integrals have been very successful in pro- 
viding predictions and cross-checks for cutting-edge analytical calculations, e.g. various 
types of two-loop box integrals or three-loop vertex functions. A restriction of the method 
for multi-loop integrals presented here is given by the fact the numerical evaluation is 
straightforward only for Euclidean points, where all kinematic invariants are negative. 
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For one-scale problems, like massless two-point or three-point functions, this is not a 
restriction at all, but if more than three external legs or/and masses are present, there 
will be branch cuts and thresholds which hinder a straightforward numerical evaluation. 
Solutions to this problem already have been suggested [32,41,42,58,59,62,64] and are 
subject to current research. 

Although the algorithm is valid to all orders in principle, there are certainly lim- 
itations from CPU-time and memory once a certain degree of complexity is reached. 
It is not possible to make a general statement about where exactly the limit is, as it 
depends not only on the computing resources but also on the way the algorithm is pro- 
grammed. Further, the number of loops and scales is not the only measure of complexity. 
Non-planar diagrams in general lead to more complicated expressions, often containing 
spurious singularities with worse than logarithmic behaviour at intermediate stages. 

In combination with its application to infrared singular real radiation, the method 
of sector decomposition has proven very useful recently to obtain differential results 
for full processes at NNLO [24,26,28,29,31,38,41]. The main advantages compared to 
methods based on analytic subtraction are the following: the subtraction procedure lends 
itself to automation and, due to the "local" nature of the subtraction terms, leads to 
expressions with good numerical behaviour. There is no need for an analytic integration 
of subtraction terms over the singular phase space regions, which is a big advantage in 
the presence of massive particles. The main drawback of the method consists in the fact 
that it leads to very large expressions for complex processes, as the number of terms is 
increased in each decomposition step. In particular for processes involving many massless 
particles, necessitating a large number of subtraction terms, like e.g. e + e~ — > 3jets or 
pp — > 2jets at NNLO, the size of the expressions produced by sector decomposition 
reaches a limit where differential results with sufficient numerical precision cannot be 
obtained within reasonable CPU times. Fortunately, most processes relevant for high 
precision phenomenology involve both massive and massless particles, where the method 
of sector decomposition has an enormous potential, not suffering from the limitations 
imposed by analytic integrability. 

Part of the problem with intractably large expressions is related to the fact that the 
algorithm, in its fully automated form, makes a decomposition already if the necessary 
condition to produce a singularity in 1/e is fulfilled (e.g. vanishing of the function T 
in the case of loop integrals). However, this is not always sufficient to produce a sin- 
gularity. Knowledge about the physical singularity structure (i.e. the soft and collinear 
limits) and inspection by eye of certain terms can certainly help to prevent unnecessary 
decompositions, but the applicability of such criteria is rather limited for complicated 
expressions as occurring e.g. in NNLO matrix elements, where a fully automated treat- 
ment is mandatory. Therefore, in order to minimise the number of produced terms, it 
would be useful to have an algorithm which finds the minimal number of decompositions 
necessary to extract the singularities. In the case of scalar loop integrals, this is basically 
a mathematical problem. If full processes are considered, a solution depends crucially 
on the way the numerator functions are treated. In any case this issue deserves further 
study. 

Finally, it is clear that the key to an optimal solution often consists in combining 
several methods in a clever way. The universality of sector decomposition, as a gen- 
eral method to isolate singularities from parameter integrals, suggests that it is a good 
candidate for such combined approaches. 
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